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1. INTRCDOCTTON 


Significant advances in present state-of-the-art analytical and 
experimental techniques are required to improve both the durability and 
reliability of hot section Earth-to-orbit engine components. These 
components, directly in the hot gas flow path, are subjected to severe 
thermal and mechanical loadings which can lead to creep-enhanced distortion 
and low-cycle fatigue. As operating temperatures are pushed higher and as 
more demanding applications, such as the space shuttle main engine and 
future Earth-to-orbit propulsion systems, are attempted, the environment in 
which the hot section components must survive becomes even more severe, and 
the effects of interaction between these components and the hot gas becomes 
significant. Consequently, the development of increasingly durable, 
structurally-efficient designs will benefit from the use of enhanced 
analytical techniques capable of producing reliable stress, deformation, 
and temperature profiles for the combined fluid-structure problem. An 
analysis of these type of problems requires the combined use of solid 
mechanics, fluid mechanics and heat transfer theories. 

This report details progress made, during the period November 1986 - 
November 1987 in a three-year program commencing in March 1986, toward the 
development of a boundary element formulation specifically designed for the 
study of hot fluid-structure interaction in Space Shuttle Main Engine 
(SSME) hot section components. The primary thrust of the program to date 
has been directed quite naturally toward the examination of fluid flow, 
since boundary element methods for fluids are at a much less developed 
state. During the first year, work focused on the completion of a 
comprehensive literature review of integral methods in fluids, the 
development of integral formulations for both the solid and fluid, and some 
preliminary infrastructural enhancements to a boundary element code to 
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permit incorporation of the fluid-structure problem. However, in the 
second year, emphasis shifted to the implementation and validation phases. 

In particular, during the past year, boundary element formulations 
were implemented in two-dimensions for both the solid and the fluid. The 
solid was modeled as an uncoupled thermoelastic median under plane strain 
conditions, while several formulations were investigated for the fluid. 
For example, both vorticity and primative variable approaches have been 
implemented for viscous, incompressible flew. More recently, compressible 
versions have also been developed. All of the above boundary element 
implementations were incorporated in a general purpose two-dimensional 
code. Thus, problems involving intricate geometry, multiple generic 
modeling regions, and aribitrary boundary conditions are all supported. 

In the next section, a brief review of the most recent boundary 
element literature for fluid flew is presented. This is followed by the 
development of integral formulations for the solid in Section 3 and for the 
fluid in Section 4. Several numerical examples are presented at the end of 
each of those two sections. In the latter case, this is followed by a 
discussion of the results and suitability of the various algorithms. The 
remaining sections then surrmarize the progress achieved to date, outline 
the work plan for the next year, and provide a list of references. Tables 
and figures appear at the end of the corresponding section. 
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2. LITERATURE PEVTEW ON FLUIDS 

Relatively little work has been done on the numerical solution of 
fluid flow problems for intermediate and large Reynolds numbers. Much of 
that work is on hybrid formulations, with integral equations used for the 
surfaces and finite differences or finite elements for the fluid volume. 
Pure integral equation formulations have the additional complexity that 
volume integrals enter the picture. This, in retrospect, is not a serious 
problem in that volume computations are restricted to those parts of the 
problem where viscous effects (both Newtonian and non-Newtonian cases) are 
significant. Furthermore, the vast majority of existing applications are 
simply for two-dimensional test problems. The following basic formulations 
are distinguished: 

Velocity and Vorticity Fomulations 

Here pressure is replaced by vorticity, which is defined as the curl 
of the velocity vector. The kinematics of flow are described by a 
vectorial Poisson's equation. The kinetic aspects of the problem are 
described by the vorticity transport equation, which is a parabolic, time- 
dependent type of equation. Integral equation formulations for either 
vectorial elliptic or parabolic equations are well understood (Wu, 1976b; 
Banerjee and Butterfield, 1981). The coupling of these two equations, 
however, is a much more difficult proposition. For the general case of 
flow past a solid object such as an airfoil, either the velocity or the 
vorticity must be specified on its surface. In addition, radiation type of 
boundary conditions must hold for unbounded fluid domains. This 
formulation is especially well suited for two-dimensional problems, because 
the vorticity becomes a scalar. 

Starting in the early seventies, WU and his co-workers have published 
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extensively on their hybrid approach for time-dependent viscous flow (Wu 
and Thompson, 1973; Wu et al 1974; Wu, 1976a; Sankar and Wu, 1978; Wu and 
Rizk, 1978; Wu et al, 1984). This work is summarized in Wu (1982; 1984) 
and the general strategy is as follows: The integral representation of the 
kinematic part of the problem is by Biot-Savart's law that involves both 
surface and volume integrals. For the case of flow past a solid body with 
the frame of reference attached to that body, the surface integrals vanish 
and what is left is a relation between the velocity and the volume integral 
of the vorticity. Viscous effects are originally concentrated only around 
the surface of the solid body and subsequently spread to a small portion of 
the surrounding fluid domain. Thus, it is only necessary to discretize the 
boundary of the solid and the part of the volume where viscous effects are 
important in order to evaluate the volume integral. Since the velocity is 
zero at the solid surface, a relationship between vorticity at the surface 
and vorticity throughout the volume is obtained through nodal collocation 
(condition 1). 

The integral representation of the kinetic part of the problem is in 
terms of both time and space integrations. Furthermore, space integrations 
include volume terms that account for the nonlinearities in the vorticity 
transport equations. Other than recasting this volume integral into more 
convenient forms, it is not in general possible to convert it to a surface 
integral. Standard nodal collocation then gives, at time step k, a 
condition involving vorticity and its spatial derivative at the surface, 
velocity-vorticity products throughout the volume, and the initial 
conditions (condition 2). 

A typical solution scheme is incremental /iterative in character. From 
known values of the velocity and the vorticity at time step k-l and 
throughout the fluid domain, condition 2 can be used to solve for the 
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vorticity at time step k and at all nodes not on the surface. An iterative 
procedure is required due to the presence of velocity times vorticity 
terms. Next, boundary values of the vorticity at k are found by using 
condition 1. Finally, interior values of the velocity at k are computed 
from the original kinetic condition. Boundary values of the velocity are 
zero for this type of problem. Using this basic type of procedure Wu and 
his co-workers solved a wide range of problems such as flow past a finite 
flat plate under zero angle of attack, flew past circular cylinders, and 
flew past airfoils. In earlier work, the vorticity transport equation was 
solved by finite differences (Wu and Thompson, 1973), finite differences 
with segmentation of the flow field (Wu et al, 1974), finite elements (WU 
et al, 1978a), and finite Fourier series in conjunction with the integral 
representation in a conformally transformed plane (WU and Sugavanam, 1978). 

A similar approach was independently proposed by Coulmy (1976). At 
first, the panel method is used for the numerical solution of the integral 
equation for kinematics. This allcws for expressing the velocity in terms 
of the vorticity. Subsequently, the vorticity is computed in an iterative 
fashion from finite differencing of the vorticity transport equation. This 
work is reviewed in Coulrty and Luu (1984), where examples of flow past 
turbine blades are presented. An approach where the vorticity field is 
discretized by means of vortex carrying particles and velocity is related 
to vorticity via Biot-Savart integral law was recently proposed for planar 
problems by Choquin and Huberson (1986). 

Vorticity and Stream Function Fornulations 

In this approach, the continuity equation is written in terms of the 
stream function, the vorticity transport equation is retained, and velocity 
is related to the vorticity in the usual way. A hybrid approach was 
proposed in Wu and Sampath (1976), where an integral equation was used to 
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represent the continuity condition. Solution was obtained in a transformed 
domain, with finite differencing used for the vorticity transport equation 
in conjunction with analytical /numerical treatment of the integral 
equation. Onishi et al (1985) and Onishi (1986) used standard integral 
equation representations for the coupled fluid flow-heat conduction system 
of equations. Heat conduction is described by a parabolic type equation, 
which is simpler than the vorticity transport equation. The nonlinear and 
coupling terms were all expressed through volume integrals. The final 
system of integral equations in both space and time was solved by the 
boundary element method. The examples solved were two-dimensional and 
included isothermal channel flow, isothermal flow past a cylinder, and 
convection-diffusion, wind driven flow in a square cavity. This work was 
an extension of the authors' work on nonlinear heat transfer problems 
(Onishi and Kuruki, 1986). 

Primitive Variable Formulations 

Several recent efforts have focused instead on the more direct 
velocity-pressure formulation. In particular, Tosaka and Kakuda (1987) 
developed a primitive variable, incompressible, integral approach by using 
an approximate time-dependent Stokes fundamental solution. In a second 
paper, Tosaka (1987) also incorporated thermal effects. However, probably 
the most important work done on this subject is that of Piva, Graziani and 
Morino (1987) where they derived the fundamental solution for the linear 
part of the Wavier Stokes equation for incompressible flow. 

Ocmpressible Flow 

The solution of compressible flow equations is done mostly through 
finite differences or finite elements. Very little information exists on 
integral equation based solutions. El-Rafaee et al (1981), for instance, 
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extended the approach for time dependent incompressible flow of Wu and his 
co-workers described previously. This was done by adding three transport 
type equations, one for each of the following dependent variables: 
dilation, enthalpy, and density. These were coupled with the vorticity 
transport equation and an augmented integral equation representing the 
kinematic aspects of flow. A conformal transformation was used to map the 
region around the airfoil to a regular region outside a circle. An 
implicit finite difference scheme was used for all the transport equations. 
The integral equation was numerically evaluated using Fourier series. 
Problems involving flow over a plate, around a circular cylinder, and 
around a Joukowski airfoil were solved by this methodology. 

Non-Newtonian Fluids 

Solutions based on integral equation formulations for non-Newtonian 
fluids have recently appeared in the literature. In particular. Bush et al 
(1984) looked at the problem of slow steady extrusion of a jet of 
Maxwellian fluid. The standard approach of integration by parts was used 
for the boundary integral equations in terms of the velocity. The non- 
Newtonian behavior was accounted for through a volume integral of pseudo- 
body forces exactly as in elasto-plastic analysis of a solid (Banerjee and 
Butterfield, 1981). The geometry of the axisymmetric problem was 
discretized via boundary elements and the numerical solution was found to 
be in good agreement with finite element results. Coleman (1984) 
approached strong non-Newtonian incompressible flow through a complex 
variable integral equation formulation for the velocity and stream 
function. The problem of slow stick-slip flow in a container of finite 
size was investigated. 

Several, perhaps obvious, points should be noted as a conclusion to 
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this section. First, the set of differential equations governing the 
behavior of the fluid are much more complicated than those pertaining to 
the solid. Thus, it is not coincidental that applications of boundary 
integral techniques for solids has reached a much more refined level, than 
that for fluids. To the investigators' knowledge, the only integral 
formulation that has appeared for time-dependent corpressible flow is the 
work done by El Rafaee, Wu and Lekoudis (1981). Even in this work, the 
vorticity, dilatation, enthalpy, and density transport equations were 
solved by finite differences. Only the kinematics were represented in 
integral form. The problems solved were also of simple test problem 
variety. 

Also, it is evident that very little has appeared in the boundary 
element literature on coupled problems. However, the present investigators 
have developed a formulation for problems of coupled transient 
thermoelasticity, which is the subject of the following section. 


8 


3. 


INTEGRAL FORMDLATION FOR SOLIDS 


3.1 Introduction 

In the present section, a surface only time domain boundary element 
method will be described for a thermoelastic body under quasistatic 
loading. Thus, transient heat conduction is included, but inertial effects 
are ignored. Formulations have been developed for three-dimensional, two- 
dimensional and axisymmetric problems (Dargush, 1987), however, only the 2D 
plane strain case is detailed below. Separate subsections present the 
governing differential equations, the integral equations, and an overview 
of the numerical implementation. 


3.2 Governing Equations 

With the solid assumed to be a linear thermoelastic medium, the 
governing differential equations for transient thermoelasticity can be 
written: 


(X+|i) 


a 2 u.. a 2 ^ 

»X i 3Xj + 4 3Xj3Xj 


( 3X+2(i) 


ae 

ax 


o 


(3.1a) 


_ ae 
at 


= k 


a 2 e 


where 


displacement vector 
© temperature 

t time 

Lagrangian coordinate 
k thermal conductivity 

p mass density 

c specific heat at constant deformation 
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(3.1b) 
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X , 4 Lame's constants 

a coefficient of thermal expansion 

Standard indicial notation has been employed with sunmations indicated 
by repeated indices. For two-dimensional problems considered herein, the 
Latin indices i and j vary from one to two. 

Note that (3.1b) is the energy equation and that (3.1a) represents the 
momentum balance in terms of displacements and temperature. The theory 
portrayed by the above set of equations, formally labeled uncoupled 
quasistatic thermoelasticity, can be derived from thermodynamic principles. 
(See Boley and Weiner (i960) for details.) 

3.3 Integral Representations 

Utilizing equation (3.1) for the solid along with a generalized form 
of the reciprocal theorem, permits one to develop the following boundary 
integral equation: 

• • 

C po ( ^ )u 3 ( ^ t) * f tG 3a* t p (X ' t) ” F pa *U p (X.t)]dS(X) . (3.2) 

where 

a,p indioes varying from 1 to 3 

s surface of solid 

u a ,t a generalized displacement and traction 

u „ - % ei T 

*« “ Itj t 2 q] 

0 ,q temperature, heat flux 

G op' F ap generalized displacement and traction kernels (Dargush, 

1987) 

C a p constants determined by the relative smoothness of s at 5 

and, for example. 
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t ' 

" J G ap (x ' t; t’ x) fc e (x ’ t) dT 
0 

denotes a Riemann convolution integral. 

In principle, at each instant of time progressing from time zero, this 
equation can be written at every point on the boundary. The collection of 
the resulting equations could then be solved simultaneously, producing 
exact values for all the unknown boundary quantities. In reality, of 
course, discretization is needed to limit this process to a finite number 
of equations and unknowns. Techniques useful for the discretization of 
(3.2) are the subject of the following section. 


3.4 Numerical Implementation 

3.4.1 Introduction 

The boundary integral equation (3.2), developed in the last section, 
is an exact statement. No approximations have been introduced other than 
those used to formulate the boundary value problem. However, in order to 
apply (3.2) for the solution of practical engineering problems, 
approximations are required in both time and space. In this section, an 
overview of a general-purpose, state-of-the-art numerical implementation is 
presented. Many of the features and techniques to be discussed, in this 
section, were developed previously for elastostatics (e.g., Raveendra, 
1984) and elastodynamics (e.g., Ahmad, 1986), but are here adapted for 
therrooelastic analysis. 

3.4.2 Temporal Discretization 

Consider, first, the time integrals represented in (3.2) as 
convolutions. Clearly, without any loss of precision, the time interval 
from zero to t can be divided into N equal increments of duration At. 

By assuming that the primary field variables, t^ and Up, are constant 
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within each At time increment, these quantities can be brought outside of 
the time integral. That is. 


N nAt 

Gpa'tpOt.t) = } tj}<X> } G pa (X-5,t-T)dr 
n=l (n-l)At 


(3.3a) 


N nAt 

F pc* u p^ x '^ ■ 1 “!! (x> I P fo (X-5.t-t)dr . 
n=l (n-l)At 


(3.3b) 


where the superscript on the generalized tractions and displacements, 
obviously, represents the time increment number. Notice, also, that, 
within an increment, these primary field variables are now functions of 
position only. Next, since the integrands remaining in (3.3) are known in 
explicit form from the fundamental solutions, the required temporal 
integration can be performed analytically, and written as 


N+l-n 

nAt 


G p„ (x -«> * J 
< 

f Gp a (X-$,t-c)dT 

:n-l)At 

(3.4a) 

N+l-n 

nAt 


F f>. <x -5> ■ J 

| Fp a (X“5 # t"t)dT • 

(3.4b) 


(n**l) At 

These kernel functions, GjJ a (x-?) and FjJ a (x-$), are detailed in Appendix B. 
Combining (3.3) and (3.4) with (3.2) produces 


N 


N+l-n 


N+l-n 

V {, V ?) '2 J [ G pa (X-{)tJ(X) - P po (X-5)uJ(X) ] dS(X) 


n=l s 


(3.5) 


which is the boundary integral statement after the application of the 
temporal discretization. 
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3 . 4.3 Spatial Discretization 

With the use of generalized primary variables and the incorporation of 


a piecewise constant time stepping algorithm, the boundary integral 
equation (3.5) begins to show a strong resemblance to that of 
elastostatics, particularly for the initial time step (i.e., N=l). In this 
subsection, those similarities will be exploited to develop the spatial 
discretization for the coupled quasistatic problem with two-dimensional 
geometry. This approximate spatial representation will, subsequently, 
permit numerical evaluation of the surface integrals appearing in (3.5). 
Hie techniques described here, actually, originated in the finite element 
literature, but were later applied to boundary elements by Lachat and 
Watson (1976). 

The process begins by subdividing the entire surface of the body into 
individual elements of relatively simple shape. The geometry of each 
element is, then, completely defined by the coordinates of the nodal points 
and associated interpolation functions. That is. 


XU) = x i (Q = tys )x iw (3.6) 

with 

c intrinsic coordinates 

shape functions 
x iw nodal coordinates 

and where w is an integer varying from one to W, the number of geometric 
nodes in the element. Next, the same type of representation is used, 
within the element, to describe the primary variables. Thus, 


u "< ■ v C)u l (3 - 7a > 

(3 - 7b) 
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in which i/^ and t£ w are the nodal values of the generalized displacement 
and tractions, respectively, for time step n. Also, in (3.7), the integer 
« varies from one to Q, the total number of functional nodes in the 
element. From the above, note that the same number of nodes, and 
consequently shape functions, are not necessarily used to describe both the 
geometric and functional variations. Specifically, in the present work, 
the geometry is exclusively defined by quadratic shape functions. In two- 
dimensions, this requires the use of three-noded line elements. On the 
other hand, the variation of the primary quantities can be described, 
within an element, by either quadratic or linear shape functions. (The 
introduction of linear variations proves computationally advantageous in 
some instances.) 

Once this spatial discretization has been accomplished and the body 
has been subdivided into M elements, the boundary integral equation can be 
rewritten as 

N M N+i-n 

V ?> V f) - 1 1 1 j [ G.„(X(;)-{m„ < i)tJ„ 

n=i m=l SL 


where 


N+l-n 

- V x<t) ' 5)N » <l)u EL jasocun ) . 


s 



(3.8) 


In the above equation, t” B and Up B are nodal quantities which can be 
brought outside the surface integrals. Thus, 


V e)u p U) 


■ 2'2 
n=l m=l 


N+l-n 

tp„ f <V X(! >-t>V c >as<xte>) 

nl 
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( 3 . 9 ) 


- u 


P<0 




N+l-n 

(Xk)-$)N 


b> 


(c)dS(X(c)) ) 


Wie positioning of the nodal primary variables outside the integrals is, of 
course, a key step, since now the integrands contain only known functions. 
However, before discussing the techniques used to numerically evaluate 
these integrals, a brief discussion of the singularities present in the 
kernels Gp a and Fp 0 is in order. 

The fundamental solutions to the uncoupled quasistatic problem contain 
singularities when the load point and field point coincide, that is. when 
r=0. The same is true of Gp a and Fp a , since these kernels are derived 
directly from the fundamental solutions. Series expansions of terms 
present in the evolution functions can be used to deduce the level of 
singularities existing in the kernels. 

A number of observations concerning the results of these expansions 
should be mentioned. First, as would be expected, F^ has a stronger level 
of singularity than does the corresponding G*p, since an additional 
derivative is involved in obtaining F^ from G*p. Second, the coupling 
terms do not have as a high degree of singularity as do the corresponding 
non-coupling terms. Third, all of the kernel functions for the first time 
step could actually be rewritten as a sum of steady-state and transient 
components. That is. 


G ap " SS ‘ 


G aP + tr<3 ip 


F^ = ® S F + ^F* 
ap *ap aP * 


Then, the singularity is completely contained in the steady-state portion. 
Furthermore, the singularity in G^ and F^j is precisely equal to that for 
elastostatics, while the G^ and F^ e singularities are identical to those 
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for potential flow. (For two-dimensions, the subscript 0 equals three.) 
This observation is critical in the numerical integration of the kernel 
to be discussed in the next subsection. However, from a physical 
standpoint, this means simply that, at any time t, the nearer one moves 
toward the load point, the closer the quasistatic response field 
corresponds with a steady-state field. Eventually, when the sampling and 
load points coincide, the quasistatic and steady-state responses are 
indistinguishable. As a final item, after careful examination of Appendix 
B, it is evident that the steady-state components in the kernels G^p and 
F^p, with n>l, vanish. In that case, all that remains is a transient 
portion that contains no singularities. Thus, all singularities reside in 

cc gg 2 2 

the G Q p and F a p components of G ap and F a p, respectively. 

3.4.4 Numerical Integration 

Having clarified the potential singularities present in the coupled 
kernels, it is now possible to consider the evaluation of the integrals in 
equation (3.9). That is, for any element m, the integrals 

/ s G^ 1_n (X(c)-5)N (ji (c)dS(X(c)) (3.10a) 

m p 

J s Fp^ 1 2 3 " n (X(0-5)N w (0dS(X(C)) (3.10b) 

will be examined. To assist in this endeavor, the following three distinct 
categories can be identified: 

(1) The point £ does not lie on the element m 

(2) The point ? lies on the element m, but only non-singular or 
weakly singular integrals are involved 

(3) The point £ lies on the element m, and the integral is strongly 
singular. 
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In practical problems involving many elements, it is evident that most 
of the integration occurring in equation (3.9) will be of the Category (1) 
variety. In this case, the integrand is always non-singular, and standard 
Gaussian quadrature formulas can be employed. Sophisticated error control 
routines are needed, however, to minimize the computational effort for a 
certain level of accuracy. This non-singular integration is the most 
expensive part of a boundary element analysis, and, consequently, must be 
optimized to achieve an efficient solution. In the present implementation, 
error estimates, based upon the work of Stroud and Secrest (1966), are 
employed to automatically select the proper order of the quadrature rule. 
Additionally, to improve accuracy in a cost-effective manner, a graded 
subdivision of the element is incorporated, especially when \ is nearby. 
For two-dimensional problems, the integration order varies from two to 
twelve, within each of up to four element subdivisions. 

Turning next to Category (2), one finds that again Gaussian quadrature 
is applicable, however, a somewhat modified scheme must be utilized to 
evaluate the weakly singular integrals. This is accomplished in two- 
dimensional elements via suitable subsegmentation along the length of the 
element so that the product of shape function, Jacobian and kernel remains 
well behaved. 

Unfortunately, the remaining strongly singular integrals of Category 
(3) exist only in the Cauchy principal value sense and cannot, in general, 
be evaluated numerically, with sufficient precision. It should be noted 
that this apparent stumbling block is limited to the strongly singular 
portions, SS F^ and ss F eft , of the F*p kernel. The remainder of F*p, 
including tr F^ and tr F^ e » can be computed using the procedures outlined 
for Category (2). However, as will be discussed in the next subsection, 
even the Category (3) ss F i j and ss F ee kernels can be accurately determined 
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by employing an indirect 'rigid body' method originally developed by Cruse 
(1974). 

3.4.5 As sembly 

The complete discretization of the boundary integral equation, in both 
time and space, has been described, along with the techniques required for 
numerical integration of the kernels. Now, a system of algebraic equations 
can be developed to permit the approximate solution of the original 
quasistatic problem. This is accomplished by systematically writing (3.9) 
at each global boundary node. The ensuing nodal collocation process, then, 
produces a global set of equations of the form 

N 

J ( [G N+1 " n Ht n ) - [F N+1_n Hu n } ) >= {0} , (3.11) 

n=l 

where 

[qN+i-iI] unassembled matrix of size (d+l)P x (d+l)Q, with 
coefficients determined from (3.10a) 

[F N+1 ~ n ] assembled matrix of size (d+l)P x (d+l)P» with coefficients 
determined from (3.10b) and c^ a included in the diagonal 
blocks 

{t n } global generalized nodal traction vector with (d+l)Q 

components 

(u n ) global generalized nodal displacement vector with (d+l)P 
components 

10} null vector with (d+l)P components 

P total number of global functional nodes 
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M 

Q ■ 5*b 

m=l 

A m number of functional nodes in element m 

d dimensionality of the problem. 

In the above, recall that the terms generalized displacement and traction 
refer to the inclusion of the temperature and flux, respectively, as the 
(d+1) component at any point. 

Consider, now, the first time step. Thus, for N»l, equation (3.11) 
becomes 

[G 1 ] {t 1 } - [F 1 ] (u 1 ) = {0} . (3.12) 

However, at this point, the diagonal block of IF 1 ] has not been completely 
determined due to the strongly singular nature of SS F. . and SS F«„. 
Following Cruse (1974) and, later, Ahmad (1986) in elastodynamics, these 
diagonal contributions can be calculated indirectly by imposing a uniform 
'rigid body' generalized displacement field on the same body, but under 
steady-state conditions. Then, obviously, the generalized tractions must 
be zero, and 

[ SS F] {1} ■= (0) , (3.13) 

where Cl) is a vector having all (d+l)P components equal to one. Using 
(3.13), the desired diagonal blocks, ss F i; . and SS F^, can be obtained from 
the summation of the off-diagonal terms of [ SS F]. Ihe remaining transient 
portion of the diagonal block is non-singular, and hence can be evaluated 
to ary desired precision. With that step completed, (3.12) is rewritten as 

IG 1 ] (t 1 ) - [F 1 ] {u 1 } = {0} . (3.14) 
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In a well-posed problem, at time At, the set of global generalized 
nodal displacements and tractions will contain exactly (d+l)P unknown 
components. Then, as the final stage in the assembly process, equation 
(3.14) can be rearranged to form 

[A 1 ] {x 1 } - [B 1 ] Cy 1 } , (3.15) 

in which 

tx 1 ) unknown components of Cu 1 } and ft 1 } 

{y 1 } kncwn components of {u 1 } and (t 1 ) 

lA^.tB 1 ] associated coefficient matrices. 

3.4.6 S9ltfti<?n 

To obtain a solution of (3.15) for the unknown nodal quantities, a 
decomposition of matrix [A 1 ] is required. In general, [A 1 ] is a densely 
populated, unsymmetric matrix. The out-of-core solver, utilized here, was 
developed originally for elastostatics from the LINPACK software package 
(Dongarra et al, 1979) and operates on a submatrix level. Within each 
submatrix, Gaussian elimination with single pivoting reduces the block to 
upper triangular form. The final decomposed form of [A 1 ] is stored in a 
direct-access file for reuse in subsequent time steps. Backsubstitution 
then completes the determination of tx 1 }. Additional information on this 
solver is available in Banerjee et al (1985). 

After returning from the solver routines, the entire nodal response 
vectors, tu 1 } and (t 1 ), at time At are known. For solutions at later 
times, a simple marching algorithm is employed. Thus, from (3.11) with 
N=2, 

IG 2 ] (t 1 ) - IF 2 ] {u 2 } + [G 1 ] tt 2 } - IF 2 ] (u 2 ) = {0} . (3.16) 
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Assuming that the same set of nodal components are unknown as in (3.14) for 
the first time step, equation (3.16) is reformulated as 

[A 1 ] {x 2 } - IB 1 ] {y 2 } - IG 2 ] It 1 } + IF 2 ) Cu 1 } . (3.17) 

Since, at this point, the right-hand side contains only known quantities, 
(3.17) can be solved for lx 2 ). However, the decomposed form of [A 1 ] 
already exists on a direct-access file, so only the relatively inexpensive 
backsubstitution phase is required for the solution. 

The generalization of (3.17) to any time step N is simply 

N-l 

lA^tx 13 ) = IB 1 ] {y^} - ^ ( [G N+1_n ]{t n ] - [F N+1 “ n ]{u n ] ) (3.18) 

n=i 

in which the summation represents the effect of past events. By 
systematically storing all of the matrices and nodal response vectors 
computed during the marching process, surprisingly little computing time is 
required at each new time step. In fact, for any time step beyond the 
first, the only major computational task is the integration needed to form 
IG W ] and IF N ). Even this process is somewhat simplified, since now the 
kernels are non-singular. Also, as time marches on, the effect of events 
that occurred during the first time step diminishes. Consequently, the 
terms containing IG N 3 and [F N ] will eventually become insignificant 
compared to those associated with recent events. Once that point is 
reached, further integration is unnecessary, and a significant reduction in 
the ccnputing effort per time step can be achieved. 

It should be emphasized that the entire boundary element method 
developed, in this section, has involved surface quantities exclusively. A 
complete solution to the well-posed linear coupled quasistatic problem, 
with homogeneous properties, can be obtained in terms of the nodal response 
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vectors, without the need for any volume discretization. In many practical 
situations, however, additional information, such as, the temperature at 
interior locations or the stress at points on the boundary, is required. 
Ihe next subsection discusses the calculation of these quantities. 

3.4.7 interior Qua nt iti es 

Chce equation (3.18) is solved, at any time step, the complete set of 
primary nodal quantities, (u N ) and {t N ], is known. Subsequently, the 
response at points within the body can be calculated in a straightforward 
manner. For any point t, in the interior, the generalized displacement can 
be determined from (3.9) with Cp a = 6p a . That is, 

N M 

u « ( 5> - 1 ' 1 t f s G el 1_n<x < m „< OdS<X(t> > 

n=l m=l m 

- u jL / S F”o 1 ' n(X(c) " 5)N tt U)dS(X(i)) . (3.19) 

m 

New, all the nodal variables on the right-hand side are known, and, as long 
as, ? is not on the boundary, the kernel functions in (3.19) remain non- 
singular. However, when £ is on the boundary, the strong singularity in 
ss 

Fp a prohibits accurate evaluation of the generalized displacement via 
(3.19), and an alternate approach is required. The apparent dilemma is 
easily resolved by recalling that the variation of surface quantities is 
completely defined by the elemental shape functions. Thus, for boundary 
points, the desired relationship is simply 

».<«> “I <3.20) 

where N w ( c) are the shape functions for the appropriate element and 
c are the intrinsic coordinates corresponding to $ within that element. 
Obviously, from (3.20), neither integration nor the explicit contribution 
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of past events are needed to evaluate generalized boundary displacements. 

In many problems, additional quantities, such as heat flux and stress, 
are also important. The boundary integral equation for heat flux, can be 
written 

N M 

q^CS) - I ( l I tj. r s E^- n (X(0-?)N„(0dS(X(0) 
n=l m=l m 

- U L / s DSi" n (X(c)-?)N u (c)dS(X(c)) ]} (3.21) 

m 

where 


Ej w (X(t)-{) 

Dj ei (X( t )-«> 


- k 


- k 


aGjJ e (x( c >-$) 

^i 

3F^(X(C)-5) 

aq 


(3.21a) 


(3.21b) 


This is valid for interior points, whereas, when ? is on the boundary, the 
shape functions can again be used. In this latter case. 


■ ywqjct) 


(3.22a) 


0N (O 
0) 

dC 


e N 

0) 



q?(5) 


(3.22b) 


which can be solved for boundary flux. Meanwhile, interior stresses can be 
evaluated from 

N M 

«!;)<«> - 1 ( } I r s E^" n (X(c )-{)N. (t)as(X(t)) 

n=l m=l m 

- u p w f s D p ij 1 " n (X(O-5)N 0) (OdS(X(O) ]} (3.23) 
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in which 


/ 


E?. 


0^j(XU)-£) = i_2v 




ad?, B(£. 

<rr — + ) - 


n. 


*k 


ij P« 


(3.23a) 


D ^ij(X(0-?) 


2pv 

l-2v 


BF 

S ij ^ 


n 

£1 


3F ei 

(— Ei 
35^ 


+ !!ii 
H* 


) - 


P 6 ij F ?e 


(3.23b) 

Equation (3.23) is, of course, developed from (3.19). Since strong kernel 
singularities appear when (3.23) is written for boundary points, an 
alternate procedure is needed to determine surface stress. This alternate 
scheme exploits the interrelationships between generalized displacement, 
traction, and stress and is the straightforward extension of the technique 
typically used in elastostatic implementations. Specifically, the 
following can be obtained 




d! • 


°ij<5> ' ‘“k.l^l.k' 5 ” * - <> 6 ijVt>uL 


(3.24a) 


(3.24b) 



u i, j(^) 


9N w N 

TT u i<0 


(3.24C) 


in which u*^ is obviously the nodal temperatures, and, 

D ijkl = X6 ij 6 kl + 2fl8 ik 6 jl * 

Equations (3.24) form an independent set that can be solved numerically for 
<jj|j(S) and u^jU) completely in terms of known nodal quantities and 
t^ w , without the need for kernel integration nor convolution. Notice, 
however, that shape function derivatives appear in (3.24c), thus 
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constraining the representation of stress on the surface element to 
something less than full quadratic variation. 

The interior stress kernel functions, defined by (3.23), are also 
detailed in Appendix B. 

3.5 Examples of Transient Thermal Stress Analysis 

3.5.1 introduction 

The complete plane strain thermoelastic formulation detailed above has 
been implemented in a general purpose boundary element code. Included are 
facilities for multiple generic modeling regions (GMR’s), sliding or spring 
interfaces, symmetry, and arbitrary time-dependent boundary conditions. 

Several examples are presented in the following subsections to 
demonstrate the validity and applicability of this boundary-only 
formulation. 

3.5.2 Sudden Heating of an Aluminum Block 

As a first example, transient heating of an aluminum block is examined 
under plane strain conditions. Ihe block, shewn in Figure 3.1, initially 
rests in thermodynamic equilibrium at zero temperature. Then, suddenly, 
the face at Y = 1.0 in. is elevated to 100°F, while the remaining three 
faces are insulated and restrained against normal displacements. Thus, 
only axial deformation in the Y-direction is permitted. Naturally, as the 
diffusive process progresses, temperature builds along with the lateral 
stresses and c 22 . To complete the specification of the problem, the 
following standard set of material properties are used to characterize the 
aluminum: 

E - lOxlO 6 psi , v = 0.33 , 

a = 13 x10 -6 /°F , 

k = 25 in. -lb. /sec. in.°F , pc = 200 in.-lb./in. 3o F . 

£ 
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The two-dimensional boundary element idealization consists of the 
simple four element, eight node model included in Figure 3.1. A time step 
of 0.4 sec. is selected, corresponding to a non-dimensional time step of 
0.05. Additionally, a finite element analysis of this same problem was 
conducted using a modified thermal version of the computer code CRISP (Gunn 
and Britto, 1984). The finite element model is also a two-dimensional 
plane strain representation, however sixteen linear strain guadralaterals 
are placed along the diffusion length. In the FE run, a time step of 0.2 
sec. is employed. 

Temperatures, displacements, and stresses are compared in Table 3.1. 
Notice that the boundary element analysis, with only one element in the 
flow direction, produces a better time-temperature history than does a 
sixteen element FE analysis with a smaller time step. Both methods exhibit 
greatest error during the initial stages of the process. This is the 
result of the imposition of a sudden temperature change. Meanwhile, the 
comparison of the overall axial displacement indicates agreement to within 
3 % for the BE analysis and 5% for the FE run. A steady-state analysis via 
both methods produces the exact answer to three digit accuracy. The last 
comparison, in the table, involves lateral stresses at an integration point 
in the FE model. The boundary element results are quite good throughout 
the range, however, the FE stresses exhibit considerable error, 
particularly during the initial four seconds. Actually, these finite 
element stress variations are not unexpected in light of the errors present 
in the temperature and displacement response. Recall that in the standard 
finite element process, stresses are computed on the basis of numerical 
differentiation of the displacements, whereas in boundary elements, the 
stresses at interior points are obtained directly from a discretized 
version of an exact integral equation. Consequently, the BE interior 
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stress solution more nearly coincides with the actual response. 

3.5.3 Tube and Fin Heat Exchanger 

Next, a more realistic problem is examined. The thermal stress 
distribution, under transient conditions, is often required to evaluate the 
durability of proposed tube and fin heat exchanger designs. Consider a 
stainless steel tube with a wall thickness of 0.050in. brazed to a 0.020in. 
gauge fin of similar material. Figure 3.2 details the geometry. Notice 
that a fillet radius of 0.015in. is assumed between the tube and fin. 

The heat exchanger is cooled continuously by a fluid at 0°F flowing 
inside the tube. It is assumed that this cooling process is of sufficient 
duration to produce zero temperature, uniformly, throughout the tube and 
fin. Then, suddenly, at time zero the outer surfaces of the tube and fin 
are exposed to a 1000°F hot gas. The convection coefficients for the inner 
and outer surfaces are 20 and 10 in.-lb./sec.in 2o F, respectively. 
Additionally, the following material properties for the metal apply: 

E * 29x10^ psi, v * 0.30, 

a ■= 9.6xlO“ 6 /°F, 

k = 1.65 in.-lb./sec.in.°F, pc * 368 in.-lb./in. 3o F . 

e 

Differences in material behavior near the braze joint are neglected. 

For the analysis one-half of a single fin is isolated. The two- 
dimensional boundary element model is depicted in Figure 3.3. Hie model 
consists of two Generic Modeling Regions (GMR's) corresponding roughly to 
the tube plus braze fillet and the fin. The tube region contains eleven 
quadratic elements and twenty-two source points, while nine elements and 
eighteen source points comprise the fin GMR. Some additional interior 
points are included to enhance the post-processing, however the solution 
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process remains boundary-only in nature. 

The resulting temperature profiles are displayed in Figure 3.4 at 0.25 
sec., 0.50 sec., 0.75 sec., and 1.00 sec. The red portions represent 
temperatures above 750°F, while blue indicates areas below 450°F. As 
expected, the thin fin, distant from the cold fluid, heats up much more 
rapidly than the tube. The most severe thermal gradients exist near the 
braze joint. Von Mises equivalent stresses are plotted in Figure 3.5 for 
points on the inner tube surface and on the fillet radius. Interestingly, 
for this particular problem, the transients are not that severe. In fact, 
the peak transient thermal stresses exceed the steady-state values by only 
a few percent. 

3.5.4 Academic Turbine Blade 

For the final thermoelastic example, a two-dimensional slice through 
an 'academic' turbine blade is examined. The three region boundary element 
model is shewn in Figure 3.6. A total of sixty-four quadratic elements are 
employed. The blade is assumed to be manufactured from stainless steel 
with the material properties defined as follows: 

E = 29x10** psi, » = 0.30, 

a = 9.6 x10” 6 /°F, 

k *= 1.65 in. -lb. /sec. in. °F, pc = 368 in.-lb./in. 3o F. 

t 

Initially, the entire blade rests unstressed at 0°F, then beginning at time 
zero all inner and outer surfaces are heated by convection. The gas 
outside the blade is at 1000°F, while the inner gas temperature is 500°F. 
The corresponding film coefficients are 2.5 and 3.0 in.-lb./sec.in. 2o F, 
respectively, except at the very tip of the blade where an outer 
coefficient of 5.0 in. -lb. /sec. in. 2o F is assumed. 

Four temperature profiles during the initial 1.6 sec. of the startup 
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are depicted in Figure 3.7, in which red represents temperatures above 
250°F, yellow 200°F-250°F, green 150°F-200°F, and blue corresponds to 
temperatures below 150°F. Meanwhile, Figure 3.8 presents the von Mises 
equivalent stress plot at 0.40 sec. Maximum stresses, greater than 40 ksi, 
are displayed in red. 

Once again, it should be emphasized that this is a boundary-only 
solution process. No volume discretization is required, and perhaps more 
importantly, steep through-the-thickness thermal gradients can be 
accurately captured. 
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TABLE 3.1 Sudden Heating of an Aluminum Block 
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FIGURE 3.4 - TUBE AND FIN HEAT EXCHANGER - TEMPERATURE PROFILES 
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4. INTEGRAL FORMULATIONS FOR FLUIDS 


4.1 Introduction 

Next, attention turns to the hot fluid. In the following, several 
alternative integral formulations are developed for both incompressible and 
compressible flow including the effects of thermal coupling. Once again, 
subsections present the governing equations, the integral equations, and a 
description of the numerical implementation. 

4.2 Governing Differential Equations for Hot Fluid Flow 

4.2.1 Time-Deoendent Compressible Flow 

Initially, the governing equations for a general compressible, 
Newtonian fluid are presented. This set will provide the basis for the 
development of the boundary integral representation. (The derivation of 
these equations can be found in standard fluid mechanics texts. See 
Yuan (1967), for example.) 

The conservation of mass in the absence of sources and sinks in the 
medium gives the equation of continuity: 


a* + !lfV 

at 3x i 


0 . 


(4.1) 


By introducing kinematics and the constitutive law for a Newtonian 
fluid with constant coefficients of viscosity, the familiar Navier-Stokes 
equations appear: 


9v. 9 vj 

p (— + Vj ^-) = U+u) 


at 

In the above. 


jS. 

ax^axj 


a 2 v 


+ n 


8Xj0Xj 


i ap 


axi 


v. velocity vector 

p pressure 


(4.2) 
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t time 

x^ Eulerian coordinate 

P mass density 

X,ii viscosity coefficients. 

For a non-Newtonian fluid, additional terms appear in (4.2). However, 
these terms can be conveniently considered as pseudo-body forces, exactly 
as done in an elastoplastic analysis of a solid. 

Next, the balance expressed by the first law of thermodynamics in 
conjunction with Fourier's law of heat conduction gives the energy equation 
as 


pc„ ( 


80 

at 


ae 


a 2 e 


+ v i ■ k srpsr ' p 


9Vj 

9Xi 


+ y 


where, again 

© temperature 

k thermal conductivity 

and 


c v specific heat at constant volume 

Y viscous dissipation. 


(4.3) 


Note that in (4.3), the thermal conductivity has been assumed constant. 
Finally, the equation of state for an ideal fluid is introduced to relate 
temperature and pressure. That is. 


p = pH© 


(4.4) 


in which 


R gas constant. 

The equations (4. 1-4. 4) represent a coupled set of six equations with six 
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unknowns, namely v^, p, p and 6. By introducing a characteristic length, 
velocity, and density, the above equations can be rendered dimensionless. 
In this case, the familiar Reynolds number (R e ) and Prandtl number (P r ) 
appear. 

Note that these equations have some strong similarities to the 
governing equations for thermocoupled solids presented in Section 3. This 
will be exploited in the development of a fluid boundary integral 
formulation. 


4.2.2 Time-Dependent Compressible Flow - Alternative Formulation 
The governing differential equations can be recast by introducing the 
concept of vorticity and dilatation, where vorticity, ok, is defined as the 
curl of the velocity. 


(1) . = p . . , 

l e igk 


dv, 


8x j 


and dilatation, e, is the divergence of the velocity, 


(4.5) 


e 




(4.6) 


In the above, e^ k is the alternating tensor. Taking the curl of (4.5), 
the gradient of (4.6), and making use of a vector identity leads to the 
following relation between velocity, vorticity, and dilatation: 


iV 


8 V x j 


0 ( 1 ) 1 . 

e ijk 



(4.7) 


Next, equations (4.2) can be reformulated in terms of vorticity and 
dilatation. First, taking the curl of (4.2) yields. 
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(4.8) 


CU) . 


i dU) i v w i 

p ( JT * v j Jx7> - * J^T. + *i - 


o U>- 


while the divergence of (4.2) produces. 


(2£ 

at 


+ v j 
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ax 


- (X+2p) 


a 2 e 

axj0Xj 


- 8 . 2 P - 

axjdxj 


+ t 


(4.9) 


Hie functions end t> in (4.8) and (4.9), respectively, collect additional 
terms involving v ^ e and p. 

Equations (4. 7-4. 9) along with (4.1, 4.3 and 4.4) completely define 
time-dependent compressible flow. While additional equations and unknowns 
appear in this alternative formulation, there are also some advantages. 
Most importantly, the kinematics and the kinetics of the problem are 
separated, with (4.7) expressing the kinematics and (4.8) and (4.9) the 
kinetics. The kinetic equations need only be considered in regions of 
nonvanishing vorticity, for the case of (4.8), and nonzero dilatation, for 
(4.9). Since, in general, the vortical region will be confined to a small 
portion of the entire fluid domain, significant reduction in computational 
effort is possible. Additionally, it should be noted that the equations 
take on a somewhat simplified and unified form. Equation (4.7) is a vector 
form of Poisson's equation for v^» while (4.8), (4.9), and (4.3) are all of 
a similar nature, representing the transport of vorticity, dilatation, and 
temperature, respectively. 


4.2.3 Time-Dependent Inco mpressible Flow 

For incompressibility, p is constant and the continuity condition 
becomes simply 



ax. 


o. 


(4.10) 
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while the equations of motion reduce to 


avi 


av. 

> ( *T + v j a^ = * a Xj axj 


ax, 


(4.11) 


For three-dimensional problems, (4.10) and (4.11) form a system of four 
equations in the unknowns and p. The equations of energy and state are 
no longer required to determine fluid motion. However, under non- 
isothermal conditions, the fluid temperatures can be obtained from (4.3) 
after the velocities are established. The exception, to this two stage 
approach, is for buoyancy driven flow in which the body forces produced by 
temperature gradients are dominant. In this latter case, continuity 
(4.10), momentum (4.11) and energy (4.3) conservation must be satisfied 
simultaneously. 


4.2.4 Time-Dependent Incompressible Flow - Alternative Formulation 
For the case of incompressible flow, the dilatation is everywhere zero 
and the formulation of Section 4.2.2 simplifies. The kinematics reduce to 


Sx^axj " e ijk axj 


a “k 


(4.12) 


and the kinetics are completely defined by the vorticity transport equation 


3 id. 00). 0V, 

" ‘ST * v j H5T - “j 


aV. 


ax 4 a 


J Xj 


(4.13) 


Now there is a system of six equations in six unknowns, v^ and o>^. Again, 
this approach has the advantage that the kinematic and kinetic aspects of 
flow are separated. It has been used by a number of BEM researchers (Wu 
and Thompson (1973), Coulny (1976)) for a number of simple test problems in 
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two-dimensions, where the vorticity becomes a scalar. 


4.3 Integral Representations 
4.3.1 Introduction 

Although in the previous section, complete sets of governing 
differential equations for various forms of fluid flew are described, from 
the practical point of view, only the velocity, pressure and temperature 
form will be sufficiently general. During the early stages of the present 
work (1986-87), the vorticity formulation was iirplemented. It was observed 
that while this formulation has some very convenient features, 
incorporation of appropriate boundary conditions for a practical problem 
becomes a difficult task. At the later stages of the current work, it may 
be possible to incorporate these vorticity integrals within a coupled 
compressible potential flow - convective heat transfer formulation to 
provide a very cost effective method for the solution of the present 
problem. Before such an approximate method can be developed, it is 
important to examine the full scale implementation of the complete 
governing equations without any approximation. 


4.3.2 Fundamental Solution 

It is convenient to rewrite the governing differential equations (4.1- 
4.4) for the velocity - pressure - temperature formulation of compressible 
viscous flow as follows: 
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a 2 v. a 2 v. 
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(4.14a) 
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(4.14b) 


where 
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f i ” " p v at ~ pv j ax . ” p v R a Xi _Re a Xi 


av^ 


0V; 


ae ap 


avj 


* - -p r -pC V • 7r~ ~pR6 — 1 + Y 

K irv at H v v 3 a X j F a X j 


r 0v i 0v i 0v k 1 0v i * 0v i 

1 M 0x j 0x i l 3 3 xk J 0x j l 3 0x j 


p " p r + P v 


P r reference density 


P v variable density 

a^ A a P !Zi n 

at v d a X j p a Xj 


p = pR0 


One of the primary requirements of developing a boundary element 
formulation is that the fundamental solution of the governing differential 
equation (4.14) must exist. These fundamental solutions can be viewed in 
same sense as the shape functions in the finite element method. For solid 
mechanics these have been very well explored. Starting with Kelvin's 
solution (1846), investigators such as Stokes, Poisson, Boussinesq, 
Mindlin, and Nowacki have provided both static and transient solutions 
which form the basis of the boundary element formulations in solid 
mechanics. It is unfortunate that workers in fluid mechanics have not 
found any use for these fundamental solutions in the infinite space and 
therefore have not made any attempt to derive such solutions. Since the 
boundary element formulations could not be developed without these 
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solutions, a substantial amount of effort was devoted in the present work 
to successively derive more and more complete solutions of the differential 
equations (4.14). As a first approximation the compressibility terms in 
(4.14) were ignored and the complete fundamental solution for a transient 
body force and a transient heat source was derived as presented below. In 
a subsequent effort these were extended to include the effect of 
compressibility, which of course, reduces the elegance of some of these 
solutions and therefore is not reproduced here. It should be noted that 
the time dependent incompressible solution given here can also be used for 
compressible flow because the body force can be adjusted to allow for the 
effects of compressibility (i.e., terms involving p v in (4.14)). The 
solution presented below is developed for the three-dimensional space. The 
corresponding two-dimensional case can be shown to formally follow the same 
procedure. 

The governing differential equation for an incompressible, viscous 
fluid can be written as: 


(i+p.) 
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Via tiie well-known Helmholtz decomposition the vectors v^ and f ^ become 


v i = v ,i + e ijk v k.j 


(4.16a) 


f i f ,i + e ijk F k,j 


(4.16b) 


with V: f = 0 and F,* .• = 0 for definiteness. Then, equations (4.15) become 

X r 1 X 9 X 
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e (v .+e., v ,x ^(v ,+e., v ) 

u + „) — q_s&LEdl ] + r — .1 e i^- Ell 1 _ 

*■ 3x^0Xj j L 0Xj3Xj J 3x^ 

r 8 (V » l~^ilJn V in, 1 ^ -i 

p L 3t J + ^,i + e ilm F m,l 0 


(4.17a) 


a(v .H* e iljii v iii.l ) 0 

8x j 

But, since 

8(e iWWi’ . 

8 *j ■ • 


(4.17b) 


equations (4.17) reduce to 


U + „> [ jf- (v (jj ) ] * » [ 5^- <» jj) * eum 8^ <V m. jj> ] * 

- 0 [ dr «=, ♦ ««. r 1 ^ ] * r- ♦ T* - * ium 

L 3x^ 3t llm 3x^ 3t J 3x^ llm 3x^ 


v . . = 0 . 

OD 


(4.18b) 


Grouping terms, equation (4.1 8a) becomes 


3 [“ 3v T 

ax. L (x+p) v ,jj + pv ,jj "P“ p at + f J 


9V_ 


+ e 


ilm 3x^ 


7 [ pV m,jj " P jf + F m 1 = 0 


(4.19) 


For generality, the bracketed terms must vanish independently. Thus, 


3v 


(X+2n)Vjj - p w - p+f = o 


(4.20a) 
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(4.20b) 


However, after enforcing (4.18b), these become 
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uv — ^ 

- P TT - P+f ■ 0 


(4.2la) 


- p ^ + F - 


(4.21b) 


Next, a unit pulse force in the i-direction applied at point x Q and time t 0 
must be decomposed. That is, let 


^ * 6(x-x 0 )S(t-t 0 )e i . 

From the properties of the delta function in three dimensional space, this 
can be rewritten as 


f i 4n K t 


in which 


y i = x i' x io 


r - 


But noting 


A i. jj = *j.ij ~ e ijk e kln^m,lj 


thus. 


h(t-t) , e. 


o r 1 x , uL 1 

L ( r ,ij e ijk e klm ( r ,lj J 


(Then, from (4.16b), the decomposition of f^ can be written 


6(t-t ) 

( r\j 


(4.22a) 
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(4.22b) 


6(t ' t o>^ 


F k “ e klm 4n 


( t\l • 


Substituting (4.22) into (4.21) yields 


9v 6 <t-t 0 )e 1 , 

--P <F>,j -0 


’ o it - P 


(4.23a) 
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(4.23b) 


The fundamental point force solution is any particular solution which 
satisfies (4.23) along with the incompressibility constraint, 


>33 


= 0 . 


(4.18b) 


First, let the scalar function, v , be identically zero. That is, let 
v(x,t;x 0 t 0 ) = 0. Ihis obviously satisfies (4.18b), and permits the 

reduction of (4.23a) to simply 


P = - 




Carrying out the differentiation, this can be rewritten as 
_ _ S(t-t Je, y i 

p(x,t;x rt ,t rt ) = s — * (-^) . (4.24) 

° 0 4nr 2 r 

Thus, (4.24) satisfies (4.18b) and (4.23a). All that remains is to find a 
solution, v m (x,t), of equation (4.23b). With that in mind, a Helmholtz 
decomposition can again be applied. That is, let 


'm 


u 


,m 


^nij u ; 




(4.25) 


with U.. = o for definiteness. Substituting (4.25) into (4.23b) produces 
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8u .m 


r 1 r » m 1»1 l 

* L u ,mll +e mi j^j , ill J “ p L ST + ^nij ST~ J 


Mt - t o’ e k .1. 

+ e mlk 4n r ,1 ' 


or after some rearrangement 


B r 0u i B r 

3x m ^’11 ” P 3t J + ^ij 0^ L pU j , 11 “ p 


0U, 6(t-t 0 ) ej 

- 2 + — 1 ] = o . 


0t 


4nr 


(4.26) 

For generality, both bracketed expressions must vanish independently. 

However, since from (4.16a), V = o, it follows that u must satisfy 

m» m 


«.ii " 0 


and consequently , 


Hi = 0 

0t ° * 


Thus, let u = 0. From the second set of bracketed terms in (4.26), 


90 n 6<t-t 0 )e, 

- P XT * Jif ^ 

This is simply a vector version of the diffusion equation. After letting 
Uj = Uej* this can be rewritten in the form of the standard scalar 
diffusion equation as 


_ 3U 

c0 'ii “at - + * = 0 


(4.27) 


in which 


c = * 
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4npr 

A particular solution of (4.27) can then be obtained by utilizing a Green's 
function approach for the diffusion equation on an infinite domain with 
sources and zero initial conditions. Therefore. 


U( 


x,t;x 0 ,t 0 ) = J G(x,t;z,T)*Mz,T;x o ,t 0 )dV(z) 


(4.28) 


where 


G(x,t;z,x) = 


e -t 2 /4 
[4jrc(t—r) ] 


3/2 


(4.29) 


is the infinite space point pulse source Green's function with 


2 _ Y 
n ~ c(t-x) 


Y = lx - z I . 

Furthermore, the volume integration in (4.2 8) is conducted over the 
infinite space v m and the operator • represents a Riemann convolution 
integral defined by 

f fc r fc 

g(t)*h(t) = g(t-r)h(r)dr = g(r)h(t-T)dr . 

J 0 J 0 

Now, also define 
e - li - x o l . 

so that the source term can be rewritten as 






6 ( T-C 0 ) 
4npc * 


Thus, 
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(4.30) 
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H(t-t 0 ) e^ 2/4c(t ~ t o ) 

4jTpe [4nc(t-t 0 )] 3/2 


with a redefinition of n as 


2 = r 2 # 

11 c(t-t Q ) 

Then, the substitution of (4.30) into (4.28) yields 


U(x,t;x 0 ,t 0 ) = } 


H(t-t 0 ) 


/ 4 


v_ 4npe [4nc(t-t 0 )] 3/2 


dV(z) . 


(4.31) 


TO carry out this volume integration, a spherical coordinate system can be 
constructed with the point x Q at the origin, as shown in Figure 4.1. 

Then equation (4.31) becomes 


U(x,t;x Q ,t 0 ) 
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Next, using the substitution 


S 4 


the integral over t\ becomes 
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New, make use of the substitutions 
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= [4c(t-t )J 1 / 2 erfc ( - 


2 Uc(t-t 0 ^ J 
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Thus, 


U(x, t; x 0 t 0 ) = 


H(t-t ) r -r i 

[ erfc( jjr) -erfc( rjr) J , 

[4c(t-t o )] 1/2 [4c(t-t 0 )] 1/2 

( 4 . 32 ) 


8npr 


but, since 


erfc(z) = 1 - erf(z) 


erf (-z) = - erf(z) 


equation (4.32) simplifies to 

H(t-t 0 ) r 

U ( x.t;x 0 ,t 0 ) - 4n P r' erf ( 


-) . 


(4.33) 


^4cTt-t“) 


With the function U now established, the fundamental solution can be 
reconstructed. From (4.16a) and (4.25), 

Vi = v,i + e ijk u <kj + e ijk e klm U mflj • 

However, since both v and u are identically zero, this reduces to 


v i e ijk e klm U m,lj ' 


which can be rewritten as 
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or, for a unit pulse point force in the j -direction. 


v i<*- t! VV - (0 aj - 8 ij U, kk >e j . 


Therefore, the desired fundamental solution becomes 


G ij< x ’ t; W ■ u ,ij - 


• _ . »<t-V y-i 

P.(x.t;x 0 ,t 0 ) - (J) 


with 


H<t-t ) 

U = erf ( - ) 


4npr 


J 4c(t-t Q ) 


(4.34a) 

(4.34b) 

(4.34c) 


in which G — is the velocity in the i-direction and Pj is the pressure, due 
to a unit pulse point force in the j-direction. 

As a first approximation, this can be combined, along with the 
standard diffusive unit pulse heat source solution to produce a boundary 
integral formulation for non-isothermal incompressible flow. In the 
compressible case, a reference density is introduced, as shown in (4.14). 
Furthermore, a dilatational component apears in G^j» along with thermal 
coupling terms. The investigators are currently developing the required 
fundamental solutions 'for the set of differential equations (4.14). Most 
of the algebraic work is completed. Some tidying up of the final solution 
is in progress. 


4.3.3 Reciprocal Theorem 

For any two unrelated states of body forces, surface tractions, 
velocities and temperatures ( f ^ 1 ^ , t ^ 1 ^ » v 1 ^ » 9 ^ 1 ^ ) and 
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existing in the same region v bounded by the 
surface s, the following reciprocal integral identity holds: 

J [ t{ 1) *vj 2) + q (1) *0 (2) - v< 1J *t' (2) - e (1) *g (2) ]dS 

S 

* | [ f{ 1, *v{ 2) + - vp’-f : <2> - e (1> *t <2) ]av , o 

V 

(4.35) 

where is used to denote modified surface tractions due to the inclusion 
of a portion of the stress tensor within f^, a s indicated in (4.14). Note 
that in equation (4.35), once again, the * indicates a Riemann convolution 
integral in time. 

4.3.4 Boundary Integral Equation 

By introducing the states corresponding to that given by 

fundamental solution (unit body force and unit heat source) in the above 
equation, one obtains the following integral representation for the 
surface velocity: 

C P a^>V 5 ' T) = J [ Sa* t P (X ' t> ' F pa* V P (X ' t} ] dS(X) 

s 

+ J [ Sa* f p (Z ’ t) J 37 * 25 (4.36) 

V 

where in two-dimensions 
v p = Ivj v 2 G] T 

, 9 f t m 

t 2 q] T 

r 9 » t m 

As in thermoelasticity, the Greek subscripts range from one to d+1. where d 
is the dimensionality of the problem. Next, perform integration by parts 


5 6 


on the third term in f^ to rewrite the boundary integral in terms of total 
tractions The result can be written as 


• • 

C ( « ( «>V 5 - T) * 1 1 G p«‘ t p (x ’ t) ‘ F p« ,v p <x - t) l 35 ®’ 


3G. . 


* vvj • • • 

J [ ^ 8 j«*Pv<2-t> + ‘W 2 ’*’ J*' z> 


(4.37) 


V 


where 


fc p ” ^ 


V 


[flf 2 ]■ 


D° 


P v = Py* 6 


1 for j = a 
0 for j # a 
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3V^ 8v^ 

p W + pv j 3x7 


The volume integral in (4.37) can be condensed further and expressed as: 

• • 

<V ?) V 5 ' T) = I [ Gp a * t p (x ' t) - Fpa’VX' 0 1 3S(X) 

s 

+ J [ VPv (Z ' t) + G pa* f P (Z ’ t) 1 3V(Z) * 


(4.38) 


For interior quantities, velocities can be computed from (4.3 8) with Cp a = 
Sp a , and additionally. 


3v . . ‘ 

“^. (5 ' T) = J [ E P ai* t p (X ' t) - D Pai* v p (X * t} F 3 


(X) 


+ I [ P ai*Pv (Z ' fc) + E pai* f p (Z ' fc) l dV(Z) + J aiPv ( ^' Tl 


V 


where is a discontinuity term arising out of the treatment of the 
improper volume integral involving P^. Meanwhile, the interior velocity 
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and temperature rates can be computed from: 


ST (5 ' T) = J t Gp 0 *t p (X.t) - Fp a *v p (X,t) ]dS(X) + j [ Gp a *fp(Z,t) ]dV(Z) . 

S V 

4.4 Numerical Implementation 

4.4.1 Surface and Volume Integrations 

The numerical treatment of the equations in fluid dynamics coupled 
with heat transfer follows very closely that described in Section 3 for 
transient thermal stress analysis. Interestingly, all of the singularities 
of the functions G and F on the boundaries, as well as, in the interior are 
identical. However, now the volume must be subdivided into cells. The 
geometry of each cell is again defined fcy nodal points and quadratic shape 
functions. In two-dimensions, six and eight-noded cells are available. 
The quadratic geometric variation permits representation of intricate 
shapes with a minimal number of cells. Meanwhile, either a linear or 
quadratic variation can be employed for the functional representation. 
Formally, then, for any cell, 

ZU> - 2.(5) - V S)Z iW 
£„( 5) = M„C 5)fJ„ 

where 

c intrinsic coordinates 

M w .M w shape functions 
2 • nodal coordinates of cell 

^ao) nodal value of body force . 
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In the above, the index w varies from one to the number of geometric nodes, 
while <o ranges from one to the number of functional nodes in the cell under 
consideration. With this spatial discretization in mind, the integral 
equations can now be rewritten. For the generalized velocity this becomes 



f s G^ 1 ' n (X(C)-5)N a) (C)dS(^)) 


/ S F^ 1 - n ( x( t )-«)N„ ( c )dS<X( c )) 1 
m 

+ / I fE„ At c5! 1 ” n (Z(c)-5)M i| (c)dV(Z(C)) 1 ) 
m=l p<0 V m Pa 


where 

L 

v = 5 v m 

m=l m 

and L is the number of volume cells. A similar expression results for the 
strain rates. The integration techniques required for the evaluation of 
the surface integrals have been discussed in Section 3.4.4. For volume 
integration, as in the surface integration techniques, subsegmentation and 
variable quadrature order are used to control error with separate schemes 
employed for singular and non-singular cases. Specific details on this 
volume integration algorithm can be found in Mustoe (1984), Raveendra 
(1984), and Banerjee and Raveendra (1986). 

The explicit form of the kernel functions Gp Q , Fp Q , and Dp ai have 
been developed for the two-dimensional case and are presented in Appendix 
C. As would be expected, these functions reduce to within a spatially 
independent constant of their steady-state counterparts when t — > ®. 
Furthermore, in a manner identical to that described for the thermoelastic 
solid, these kernels decompose into a singular steady-state portion plus a 
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non-singular transient component. 


4.4.2 The Treatment of Continuity Equation 

The continuity equation for a compressible flew takes the role of the 
elastoplastic constitutive equation in solid mechanics. Recalling the 
continuity equations 


p 


p r + Pv 


(4.39a) 


at 


V j 



av^ 

p -r^- - ° 
axj 


(4.39b) 


where p is the density. 

It is necessary to satisfy these equations at discrete nodal points. 
By defining a global shape function, the density variation can be 
interpolated as 


M 

p(x X ) = Y C(x 1 ,? m )©(5 m ) for 1 = 1,2,...,M (4.40) 

m=l 


where C(x^,? m ) = l — - — where r = (y^yj) 1 ^ 2 and y. = 

r max 1 ill 

r max = rnax i mum distance between any two points of the region. 

In order to satisfy (4.39), the values of p (or p v ), dp/dxy Vj and avj/axj 
must be calculated. Of these, v^ and av^/axj are obtained directly from 
the boundary or interior integral representations. The value of p or p v is 
obtained from (4.40) and its gradient is calculated via 

M 

f^x 1 ) = ^ Dj(x 1 ,? rn )0(? m ) (4.41) 

m=l 

where 


60 








By utilizing the known values at time t in the continuity equation, the 
unknown density at t+At can be obtained assuming a Crank-Nicolson type y- 
algorithm. Thus: 


[L lm J <V = 


{ 1 ^} 


(4.42) 


where 


L lm " 


[ It ♦ * t+At 


9v • 


(X 1 ) ^ * [ r^W) W 


im 


*1 - [ it ♦ <n»‘ £ O' 1 ’ fp * [ (r-i-V 1 ’ 


Once (0 m ) is determined, the required quantities t+At p(x^) and t+At |£-(x^) 

j 

can be computed from (4.40) and (4.41). 


4.4.3 Algorithm for Incremental Iterative Solution 
The final system of algebraic equations resulting from the previous 
integral representations can be expressed as 


A b x 

- (Pf = Bhy 

(4.43a) 

V = 

A v x + B v y + G v f 

(4.43b) 

• 

V « 

A v x + B v y + G v f 

(4.43c) 

c = 

A e x + B c y + G 8 f 

(4.43d) 


where 

x,y are the known and unknown boundary quantities 
v.v are the velocity and acceleration vectors 
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e are velocity gradients 


An iterative algorithm similar to the initial stress method (Banerjee and 
Butterfield, 1981) can then be developed as follows: 

1. Increment boundary conditions 

2. Assume f = 0 

3. Calculate the boundary and interior solutions x,v,v and e 

4. Determine f, K and p v at this time increment 

5. Calculate the boundary and interior solutions again 

6. If the solution is not significantly different from (3), go to 
(l); if the solution is different, then go to (4). 

4.4.4 Algorithm for Direct Solution 

Just as the finite element method can be used to obtain the solution 
of a nonlinear system either iteratively or directly using incrementally 
updated system equations, analogous techniques can be developed for the 
boundary element method. The direct BEM algorithm was developed recently 
for elastoplasticity by Banerjee and Raveendra (1987) and Henry and 
Banerjee (1988). For the computational fluid dynamics problem the same 
underlying principle can also be used to develop a direct algorithm. 
Equations (4.43) can be rewritten for this purpose as: 

A b x + Bky + G b f = 0 (4.44a) 

e - A e x - B e y + G e f = 0 (4.44b) 

By pre multiplying (4.44b) by v T one obtains 

V T e - v T A e x - v T B e y + V T G e f = 0 . (4.45) 

Noting that v T e = f, (4.45) and (4.44a) combine to produce 
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(4.46) 
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Equation (4.46) can be solved directly if v, p y or p can be estimated for 
the next increment. This does, of course require at least one iteration. 

If substructuring is necessary, the above equation can be expressed 
for a problem involving n subregions as: 
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v T A, 


( I- V^G j 
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v T A 2 d-v T G f ) 
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0 A^ 0 
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V -°n 
v\ (i-v T Gf ) 



This block banded set of equations must be formulated and decomposed 
on each iteration. 


4.5 Examples of Viscous Plow Analysis 
4.5.1 Introduction 

The various fluid formulations, described in the foregoing, are in 
various stages of development. For incompressible flow, both iterative and 
direct algorithms have been implemented, although validation has been 
completed only under steady-state conditions. On the other hand, the 
compressible flow formulation has again been implemented, but serious 
testing of the algorithm awaits completion of the new time-dependent 
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compressible kernels mentioned in Section 4.3.2. 

In the following subsections, several examples are presented to 
demonstrate the validity of the present boundary element formulations. 

4.5.2 Driven Cavity 

The two-dimensional driven cavity has become the standard test problem 
for incompressible computational fluid dynamics codes. In a way, this is 
unfortunate because of the ambiguities in the specification of the boundary 
conditions. However, numerous results are available for comparison 
purposes. 

The incompressible fluid of uniform viscosity is confined within a 
unit square region. The fluid velocities on the left, right and bottom 
sides are fixed at zero, while a uniform non-zero velocity is specified in 
the x-direction along the top edge. Normal (y-direction) velocity is also 
zero on the top surface. Thus, in the top corners, the x-velocity is not 
clearly defined. To alleviate this difficulty in the present analysis, 
this velocity component is smoothed out as indicated in Figure 4.2. 
Additionally, steady-state conditions have been assumed and the forcing 
velocity is slowly incremented to achieve the desired Reynolds number. 

Several levels of mesh refinement, as shown in Figure 4.3, were 
examined to determine the effect on accuracy and convergence. Spatial 
plots of the resulting velocity vectors for the twenty-four cell model are 
presented in Figure 4.4 for Reynolds numbers (Re) of 0, 100 and 250. 
Meanwhile, Figure 4.5 contains the velocity distribution along the vertical 
cavity centerline. All of these results correlate closely with those 
obtained by Burggraf (1966) in his classic paper on steady separated flows. 
However, as Re increases, the convective terms begin to dominate the 
viscous ones, and then problems with numerical instability arise. In the 
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current boundary element analysis, results diverge as Re reaches a certain 
value dependent, initially, upon mesh refinement. Beyond a certain level 
of mesh refinement, though, no improvement was found. Table 4.1 contains 
the approximate Reynolds number at divergence for the various cell patterns 
displayed in Figure 4.3. Additionally, decreased increment size and/or 
increased number of iterations had little effect on the divergent character 
of this problem. 

The thirty-six cell model was also run as a three region problem as 
defined in Figure 4.6. This dramatically reduced computing time, since now 
integration is limited to within each individual region. At low Re, 
results were quite good, thus verifying the multi-GMR approach. 
Convergence was enhanced to the extent that solutions could be attained up 
to Re = 300. Figure 4.7 contains velocity vector plots at Re = 250, 300, 
and 325. The latter is an unconverged result. Notice the disturbances 
entering the solution near the center of the cavity. 

Next, the direct algorithm described in Section 4.4.4 was utilized for 
the same problem. In this case, computational effort increased 
dramatically, since interior forces are added to the list of primary 
unknowns. While convergence was generally achieved with only a few 
iterations in the lower Reynolds number range, the divergent character of 
the solution process remained unaltered. That is, even with a fine mesh, 
divergence occurred at approximately Re = 300. 

Interestingly, the same phenomenon occurs in both finite difference 
and finite element approaches in approximately the same Reynolds number 
range. In these methods, upwinding schemes (Allen and Southwell, 1955; 
Henrich and Zienkiewicz, 1979) and artificial dissipation (Hughes et al, 
1979) have been devised to enhance convergence. Ihis would occur almost 
automatically in a transient boundary element formulation because these 
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upwinding schemes and dissipations are part of the transient solutions. As 
an initial test of the transient formulation, the cavity was examined with 
an abrupt jump in velocity along the upper edge at time zero. Results at 
four different time steps, for R e “ 100, are shown on Figure 4.8. Notice 
that eventually the steady-state solution is achieved. 

4.5.3 Buoyancy Driven Flow 

The domain is again the unit square cavity. Velocity is zero on all 
boundaries, while the top and bottom edges are insulated. The temperature 
of the right edge is maintained at zero, but the left face temperature is 
slowly elevated under steady-state conditions. By including the body 
forces due to buoyancy, the existing temperature gradients produce free 
convection within the cavity. The resulting steady-state velocity vector 
distribution for a Grashof number of 100 is presented in Figure 4.9, while 
the corresponding temperature variation along three vertical planes is 
plotted in Figure 4.10. Again, in a manner similar to that encountered in 
the previous example, convergence difficulties arise in the steady state 
case as the temperature gradients (i.e., the Grashof number) are elevated 
beyond a critical value. A transient solution would probably not have the 
same difficulty. 

4.5.4 Unconfined Flow Around Obstacles 

Boundary element models have been constructed for flow around a 
simplified blade as shown in Figure 4.11. Preliminary results have been 
obtained. This problem will be investigated more thoroughly in the next 
few months. 

4.6 Discussion 

Obviously, based upon the results of the numerical studies presented 
in the previous subsection, significant improvements in the current 
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boundary element implementation for the fluid are still needed. First of 
all, the recent introduction of the time-dependent kernel functions 
described in Section 4.3 should help considerably. With time dependence, 
as Reynolds number increases the transition is from parabolic to hyperbolic 
behavior, rather than from elliptic to hyperbolic. From a numerical 
standpoint, this is much more forgiving. In fact, it is not surprising 
that most finite difference and finite element approaches to steady-state 
flews utilize either a real or psuedo time stepping approach. 

Secondly, the introduction of compressibility into the system should 
have a positive effect. Experience in plasticity, for example, has 
indicated that more accurate solutions are attainable for compressible 
materials. In addition, there is considerable room for improvement in the 
non-linear solution algorithms. Presently, only simple iteration has been 
investigated, however, implicit schemes, Newton-Raphson methods, and 
optimization techniques may prove beneficial. 
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FIGURE 4.8 - DRIVEN CAVITY - TRANSIENT RESPONSE 
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5. SUMMARY CP PROGRESS 

As is evident from the previous two sections, significant progress has 
been made toward the goal of developing a general purpose boundary element 
method for hot fluid-structure interaction. For the solid phase, a 
boundary-only formulation was developed and implemented for uncoupled 
transient thermoelasticity in two dimensions. The elimination of volume 
discretization not only drastically reduces required modeling effort, but 
also permits unconstrained variation of the through-the-thickness 
temperature distribution. Consequently, it is expected that this 
formulation can provide an attractive alternative to finite element methods 
for transient thermal stress analysis, particularly when accurate stresses 
are required. It should be noted that the BEM detailed in Section 3 is the 
first boundary-only implementation for this class of problems. 

Meanwhile, for the fluid, significant strides have also been made. 
Fundamental solutions have been derived for transient incompressible and 
compressible flew in the absence of the convective terms. Boundary element 
formulations have been developed as well and were described in Section 4. 
For the incompressible case, the necessary kernel functions, under 
transient and steady-state conditions, have been derived and fully 
implemented into a general purpose, multi-region boundary element code. 
Several examples have been examined in detail to study the suitability and 
convergence characteristics of the various algorithms. 

While considerable effort is still needed to produce a useful 
analytical tool for high Reynolds number hot gas flow, the development to 
date does represent a substantial advancement in the state-of-the-art. 
Table 5.1 is provided in order to highlight the status of the major 
components of the overall hot fluid-structure interaction program. 
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6. WOFKPLAN FOR THE NEXT YEAR 

The following tasks are planned, in approximate chronological order, 
for the period November 1987 - March 1989: 

1. Complete the time-dependent compressible flow kernel development 
and general-purpose implementation. 

2. Begin validation of the time-dependent formulations. After the 
initial debugging phase, download the code to the NASA Cray for 
this task. 

3. Investigate enhanced non-linear solution algorithms for Navier- 
Stokes flow. 

4. Develop and implement a thermally-sensitive inviscid, 
compressible flow formulation. 

5. Implement interface coding for fluid-solid interaction. 

6. Assess suitability for large scale development of boundary 
element methods for hot fluid-structure interaction. 
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APPENDIX B - KERNELS FOR TBERMDEEASTICm 


This appendix contains the detailed presentations of all the kernel 
functions utilized in the formulations contained in Section 3. IVo- 
dimensional (plane strain) kernels are provided, based upon continuous 
source and force fundamental solutions. For time-dependent uncoupled 
quasistatic thermoelasticity the following relationships must be used to 
determine the proper form of the functions required in the boundary element 
discretization. That is, 

GjJp(X-5) = G ap (X-S,nAt) for n=l 

gJJ p <x- 5) = G ap (X-5,nAt) - G a p(X-£, (n-DAt) for n>i , 

with similar expressions holding for all the remaining kernels. In the 
specification of these kernels below, the arguments (X-£,t) are assumed. 
The indices 

i,j,k,l vary from 1 to d 

a,p vary from l to (d+l) 

© equals d+l 

where d is the dimensionality of the problem. Additionally, 
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APPENDIX C - KERNELS FOR VISCOUS, INOOMFRBSSIBIJB FLOW 


This appendix contains details of the time-dependent incompressible 
kernels necessary for the integral formulations of Section 4. Notation is 
consistent with that defined in Appendix B. 

For the generalized velocity kernels. 
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